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We develop a robot algorithm to maximise the number of distinct states reliably extracted from 
correlator data using the variational analysis method. The robot explores the variational param- 
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1. Introduction 

The variational approach [1] is a widely used method to extract the excited state spectrum from 
matrices of coiTclation functions. Apart from the choice of operator basis, the variational method 
has a parameter space that should be expected to yield consistent results for a range of values 
(within certain limitations [2]). Exploring this parameter space potentially yields lai^ge number of 
correlation functions to fit. 

We also currently have access to accurate methods for the determination of all-to-all quark 
propagators, e.g. Refs [3, 4, 5, 6, 7], which allows easy access to a very large number of available 
basis states for the variational approach. This in turn leads to another multiplicative factor of 
con^elation functions to be fitted. Ultimately we must automate the procedure for fitting these data 
sets since, at some point, it becomes unfeasible to do it by hand. 

Here we describe both these processes and formulate a fitting procedure that attempts to com- 
bine a fitting algorithm and an exploration of the variational parameter space. It strives to deliver 
figures from correlation matrix analysis which can be used to determine consistent and robust fitted 
energies for multiple energy levels. 

1.1 Variational Analysis 

Extracting excited-states using a variational approach requires a matrix of correlation functions 

c«/3(0 = (o|o«(0oj(o)|o) 

where Oa (oc = I ■ ■ .N) form the operators with the appropriate quantum numbers of our varia- 
tional basis. We define the A'^ principal correlators Xa {t,to) as the eigenvalues of C(?o) '^^C(f)C(fo) 
(where to is the time defining the "metric"). tQ is (in theory) sufficiently large that one believes that 
the correlator at is dominated by the lightest N states and (in practice) sufficiently small that the 
inversion of C(?o) does not become numerically unstable. We can then show [1] that 

limA„(?,?o) = ^'"^'"'°'^°(l+0(^'""^°)) where AEa = minlEg - Ea\ . 

The N principal effective masses defined by in^^ {t) = ^n^TJi+Ua)) "^^^ '■^'^^ (plateau) to 
the A'^ lowest-lying stationary-state energies. Since all the time dependence of the diagonalisation 
process is contained in the principal correlators we can perform the diagonalisation at a single 
"optimising" time, t = topt, and the eigenvectors determined at this timeslice can be used to diag- 
onalise the correlation matrix at all times. We can also then extract the effective masses by fitting 
the plateau regions of the N diagonal entries of the projected correlation matrix. 

This method relies on constructing a basis of operators Oa that provides a good description of 
the states of interest. With the advent of all-to-all propagator techniques we can easily construct 
a fuller variational basis from a combination of the multiple operators for each Oh representation 
state (for example) and from the application of different smearing operators to the quark fields. 

1.2 All-to-all Propagators 

Many algorithms have been devised to approximate all-to-all propagators, for example see 
Refs. [3, 4]. Here we use an exact implementation involving a hybrid method that combines an 
eigenvector decomposition with a variance-reduced stochastic estimator [5]. 
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To construct the all-to-all propagator, the lowest N^y eigenmodes (eigenvectors v^'^ with eigen- 
values A,) of the hermitian Dirac matrix Q = y^M are first computed, and a truncated spectral 
decomposition of the propagator is then given by 



(1.1) 



If is equal to the dimension of the matrix, then Qq = Q ^ otherwise the propagator can be 
expressed as 

Q-' = Qa + Qu (1-2) 

and the truncation in the eigenvector representation can be corrected by estimating Q\ stochasti- 
cally. 

We choose to "dilute" the noise vectors, Tj , used in the stochastic estimation, which results in 
rapid variance reduction. In this context, dilution means creating a set of noise vectors by applying 
a set of masks to a single noise source. These masks might for example select a particular time- 
slice of the vector (referred to as time dilution), thus returning Ni noise vectors from each single 
noise source. 

The method can be implemented efficiently in software by use of a "hybrid list" method. Two 
lists of A^HL = A'ev +A'^dii vcctors u and w are written. 



w 



(0 



■^1 ' ' 



('■) = Jv(l) ... v(^ev) „,(1) 



(1.3) 



where Mi//^^'^ = 75 (/ — Pq) and Pq = ^f" v^'^ ®v^^^ . The inverse of M can then be written as a 
single sum over the pair of lists, M^' = ® h'^''^)75 . In all data contained in this paper 20 

eigenmodes and dilution in time and colour space is used. 
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Figure 1: Correlation (left) and effective mass (right) plots for the pion ground state. 

A sample correlation function, C{t), for a pion on 96 configurations of x 80 dynamical 
anisotropic lattices can be seen in Fig. 1 (these lattices are used in all results shown and the lattice 
parameters are detailed in Ref. [8]). What is interesting to note however is that the use of time 
dilution has a dramatic consequences for the effective mass plot shown on the right of this figure 
(defined here for a meson with periodic boundary conditions by / og ( ^'^'^ ^■yc(t)'^ ^ ^ ) ■ "^^^^ 
cause we introduce independent stochastic noise at each timeslice. It is the correlation function that 
is fitted however and fitting between timeslices 21 to 37 yields the fit (with en^ors) in the effective 
mass plot on the left of Fig. 2 . However, as can be seen on the right of Fig. 2 , if we allow the 
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Figure 2: Enlarged version (left) of the effective mass plot for the pion ground state using the standard 
definition. On the right is an identical plot using the a 3-timeslice extended definition. The fit between 
timesUces 21 to 37 is shown by the black line with errors given by the dashed lines. 



effective mass plot to have more extended information by using next-to-nearest (or next-to-next- 
to-nearest) neighbours we can clean up the effective mass plot dramatically and our fit starts to 
look less outrageous. The reason for this is that the fit has information across an extended temporal 
region due to its fit window and can see through the temporal stochastic noise introduced by the 
stochastic estimator. 



2. A Fitting Robot 

The variational approach has a parameter space that can be explored by variation of the "met- 
ric" and "optimising" timeslices. Also the use of all-to-all propagators can dramatically extend the 
size of the operator basis. Exploring these possibilities can lead to the creation of an extremely large 
number of correlation functions for fitting. At some point it becomes unfeasible to consider fitting 
these resultant correlation functions by hand. In the particular case detailed above, the legitimacy 
of using the human eye and an effective mass plot to ascertain the accuracy of a fit that numer- 
ically gives a reasonable reduced value is immediately brought into question. It is for these 
reasons that we have investigated the creation of a fitting algorithm to fit correlation functions. The 
algorithm is based around 3 axioms, respected in the following order, 

• A reduced value below a given minimum signifies an acceptable fit, 

• The largest possible temporal window is to be chosen, 

• Plateaus beginning at early timeslices are preferred. 

The first axiom is a standard fitting criterion applied to this type of data (here we use a value 
of 1 . 1 as the cut-off). The second axiom has a caveat in that only correlator data having a fractional 
error less than a chosen value (in the cases shown here 50%) can be included in the fit. We do this 
because we must identify the limitations of our data; including terms that do not constrain the fit 
only serve to artificially lower the reduced value. This then chooses the maximum time value 
allowed in any fit region. The third axiom arises from the fact that the lower correlator values 
heavily constrain the fit (since the signal to noise ratio is decaying). If a large fit window exists that 
includes these points then they should be included to gain the full accuracy the data allows. Another 
reason for the third axiom is that the projected correlation matrix is not, in practice, completely 
diagonalised and it is possible that we can have contamination from lower-lying states. This can 
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lead to a decay of the correlator into the lower state at large temporal values, dramatically effecting 
the accuracy of a fit. 

In practice, the algorithm is implemented by first cutting the data at the appropriate temporal 
value (the last one with acceptable fractional error). Following our second axiom, we chose the 
largest possible temporal window and attempt a fit. Unless we have managed to exactly decompose 
the correlation matrix, such a fit will fail due to contamination from higher lying states. These states 
have greatest effect at low temporal values since they decay exponentially. Since such values are 
the most accurately determined (due to their signal to noise ratio), the inclusion of these low-lying 
points will always enlarge the reduced X'^ value because they are so heavily constraining. It should 
be noted that this in effect also allows us to chose the lowest time value allowed in the fit region to be 
as low as possible (but avoiding the contact term): any fit including early contaminated timeslices 
is destined to be rejected by virtue of the first axiom. Next we decrease the allowed fit window size 
and attempt all possible fit windows in the allowed region. We chose the ordering of these attempts 
such that the starting timeslice, f„„„, of our fit window moves forward in time, i.e., we begin with 
tmin at its lowest allowed value and then increase it until we have exhausted all possibilities. This 
then satisfies our third axiom. We iterate this procedure until we achieve an acceptable reduced X'^ 
value thus satisfying our first and most important axiom. 



2.1 Results 

In the results shown below, the operators employed are a combination of the point and extended 
operators detailed in Ref. [8] and 4 different quark smearing levels (with a stout-link smeared gauge 
field). We use a quark mass in the region of the strange quark mass for all data. 
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Figure 3: Plots of the fitted effective masses and eigenvalue masses as determined for a range of metric 
timeslices for the pion. The subdivisions correspond to different optimising timeslices. 



On the left of Fig. 3 we show the fitted effective masses of the pion, determined using the 
fitting robot from the projected coiTclation functions, for a single metric timeslice, 2, and a range of 
optimising timeslices. Along the X-axis we have the metric timeslices, ?,„e,, with the sub-divisions 
between the points con^esponding to different optimising timeslices, the first point being (f„,^, + 1)> 
the second + 2), etc. We can directly observe the contamination of the eigenvectors by higher- 
lying states in the excited channel for low optimising timeslices. It is the plateauing from below of 
this effect as the optimising timeslice is increased that is interesting. For comparison we have also 
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included the mass with error determined from the eigenvalues of the variational analysis. We can 
see that the eigenvalues are indeed approaching the same plateau regions. 

On the right of Fig. 3 we compress the data of the left hand side and include further data from 
subsequent metric timeslices. For the ground state the fitted mass is consistent in all cases. The 
plateau effect in the excited channel is apparent at low metric timeslices but almost disappears from 
metric timeslice 7 onwards. We observe a consistency of these plateau regions across all metric 
timeslices. We also note that the eigenvalue effective masses begin to fall directly on the fitted 
regions as the metric timeslice increases. While the errors on these eigenvalues increase, the errors 
determined from the projected data remain unaffected. 







++ 



■>++ 









Fitted 


State ^ — 1 








Fitted 


State 1 -X- -■ . 








Eigenvalue 


State !- -«- -: 


i 






Eigenvalue 


State 1 - H -. . 








Eigenvalue 


State2 i -■ -; 












ISs S 


'i - 








%A « 








f - 






















2 3 


4 


5 6 


7 8 


9 1 






Metric Timeslice 







0.8 

0.7 ■ 

J- o.e ■ 

I 0.5 ; 



Fitted - State ^ 

Fitted - State 1 ■- 

Fitted - State 2 ■■ 

Eigenvalue - State ; 

Eigenvalue - State 1 ■ 

Eigenvalue - State 2 ; 



ft 



Metric Timeslice 



Figure 4: As before, for the 0++ (left) and 2++ (right) states. 



These operators above were just different smearing combinations of the standard 75 pion op- 
erator since this is known to have a high signal-to-noise ratio. We have also carried out similar 
analysis for a range of states and having included extended operators. Here we show some of the 
more interesting results. In Fig. 4 we can see similar behaviour to the pion for the 0++ isovector 
state and the 2++ P-wave state. Only slight variations in the overall behaviour occur in both cases. 
For the we see that the effective masses determined from the eigenvalues drops in a few cases 
at timeslices 5 and 6. At metric 6 this is accompanied by a disturbing drop in the determined ground 
state energy. However, since the overall trend is obvious, we can still accurately chose fit values 
for the ground and first-excited states. In the case of the P-wave we observe slightly more eiTatic 
behaviour of the fit to the projected data, coupled again with eiTatic behaviour of the eigenvalues. 
In this particular case, three operators were used in the con^elation matrix but after the first few 
very early successes the robot was unable to find suitable fits of the projected data. Resulting errors 
from poor diagonalisations may have caused the effects that we observe in this case. 

Lastly, in Fig. 5 we also include results for the p meson where we have large data sets for 3 
operators. On the left we have the fitted mass values; on the right we have the same data but with 
the eigenvalue masses overlaid. In the cases where we were unsuccessful in the diagonalisation 
of the 3x3 correlation matrix we have reduced it to a 2 operator diagonalisation (for example 
all of metric 9). On the left, again we see some erratic behaviour in plenty of cases but at mid 
range metrics and higher optimising timeslices we see quite systematic behaviour and accurate 
determination of 3 states. Interestingly, on the right graph, it appears that the lowest eigenvalue is 
plateauing at the first-excited state fitted mass (even in the case of metric 9), and that none of the 
eigenvalues approach the ground state value that is determined so consistently across all cases. 
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Figure 5: As before, the figure on the left is plotted without eigenvalue masses. 



3. Conclusions 

The fitting robot functions consistently in our analysis. With any such attempt however, it 
is inevitable that numerical situations will arise that circumvent all axioms of the robot due to a 
spuriously low reduced value. It is for this reason that the method couples so well with the 
variational approach. This allows for multiple fittings of data derived from an identical source and 
we can then observe plots such as those in Fig. 4 and use the human eye as the ultimate test to 
determine where the fit truly lies. Improvements to the robot should include a consistency check 
of the chosen window - to identify incorrect fits a bootstrap of the chosen data points should be 
performed to search for excessive chi-squared and fitted mass variations. To assist with this, a fit to 
the principal correlators for each metric timeslice should also be included on the derived figures. 

This work has yielded some very interesting insights into the variational approach. These re- 
sults indicate that the robustness of the choice of variational parameters should always be checked. 
It also complements Refs. [2, 9] in their discussions and conclusions with respect to the selection 
of the metric timeslice. We believe that with the suggested improvements to the fitting robot and a 
larger and more probing selection of the operator basis that such an approach can yield consistent 
and robust results from the analysis of correlation matrices. 
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